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Abstract:  Autonomous  precision  placement  of  parafoils  is  challenging  because  of  their  limited 
control  authority  and  sensitivity  to  winds.  In  particular,  when  wind  speed  is  near  the  airspeed, 
guidance  is  further  complicated  by  the  parafoils  inability  to  penetrate  the  wind.  This  article  speci¬ 
fically  addresses  the  terminal  phase  and  develops  an  approach  for  generating  optimal  trajectories 
in  real-time  based  on  the  inverse  dynamics  in  the  virtual  domain.  The  method  results  in  efficient 
solution  of  a  two-point  boundary-value  problem  using  only  a  single  optimization  parameter 
allowing  the  trajectory  to  be  generated  at  a  high  rate,  mitigating  effects  of  the  unknown  winds. 
It  is  shown  through  simulation  and  experimental  results  that  the  proposed  algorithm  works  well 
even  in  strong  winds  and  is  robust  to  sensor  errors  and  wind  uncertainty. 

Keywords:  optimal  control,  parafoils,  predictive  control,  trajectory  optimization 


1  INTRODUCTION 

Precision  payload  delivery  using  guided  parafoils  has 
extended  military  re-supply  capability  by  providing 
accurate  and  rapid  response  at  substantial  offsets  with 
a  minimal  risk  to  the  cargo  delivery  aircraft.  Signifi¬ 
cant  efforts  have  been  underway  to  improve  precision 
airdrop  with  the  US  Department  of  Defense’s  Joint  Pre¬ 
cision  Airdrop  System  (JPADS)  programmes  [1]  being 
an  example.  The  JPADS  self-guided  programmes  are 
separated  into  five  classes  according  to  size:  micro 
light  weight  (MLW)  4-70  kg,  ultra  light  weight  (ULW), 
100-320  kg,  extra  light  (XL)  320-1000  kg,  light  (L)  2200- 
4500  kg,  and  medium  (M)  6800-14  000  kg  [2].  Most 
existing  efforts  have  focused  on  larger  systems  above 
200  kg  including;  Onyx  from  Atair  Aerospace,  Panther 
from  Pioneer  Aerospace,  Screamer  from  Strong  Enter¬ 
prises,  and  Spades  from  Dutch  Space  all  demonstrated 
at  the  Precision  Airdrop  Technology  Conference  and 
Demonstration  in  2007  [3].  Recently,  there  has  been 
a  new  focus  on  smaller  systems  with  the  start  of  the 
Joint  Medical  Distance  Support  &  Evacuation  ( JMDSE) 
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programme  in  2009.  The  JMDSE  programme  is  focused 
on  small  medical  equipment  delivery  in  the  MLW  and 
ULW  classes  and  the  integration  into  small  unmanned 
autonomous  systems. 

Autonomous  parafoils  for  precision  delivery  are 
unique  aerospace  systems  with  challenges  includ¬ 
ing  limited  control  authority  and  sensitivity  to  winds. 
Parafoils  typically  only  have  lateral  control  using 
asymmetric  brake  deflection,  with  longitudinal  con¬ 
trol  using  symmetric  deflection  being  very  limited.  The 
absence  of  thrust,  combined  with  low  airspeed,  makes 
wind  a  significant  contributor  to  a  parafoil’s  trajectory. 
Airspeed  for  parafoils  can  range  from  6  m/s  for  very 
small  lightly  loaded  systems  to  over  20  m/s  for  larger 
systems.  This  results  in  the  common  occurrence  for 
parafoils  where  the  wind  speed  may  be  near  or  exceed 
the  airspeed.  In  the  case  where  wind  speed  exceeds 
airspeed,  the  parafoil  cannot  advance  towards  the  tar¬ 
get  if  the  target  is  upwind.  This  becomes  increasingly 
likely  for  parafoils  in  the  MLW  class. 

Guidance  strategies  for  parafoils  vary  significantly 
over  the  spectrum  of  systems;  however,  there  are  some 
common  aspects  seen  in  all  systems.  One  common 
aspect  is  the  parafoil  is  released  upwind  of  the  target 
at  a  sufficient  altitude  to  ensure  that  the  target  can 
be  reached.  In  addition,  all  strategies  have  some  com¬ 
bination  of  homing  and  energy  management  phases 
followed  by  terminal  guidance.  Homing  refers  to  the 
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process  of  bringing  the  system  from  one  region  to 
another  region  while  energy  management  refers  to 
loitering  in  a  pattern  allowing  the  system  to  reach 
an  altitude  where  the  parafoil  can  proceed  to  impact 
the  target.  Assuming  homing  and  energy  management 
stages  bring  the  parafoil  to  a  location  where  arriving  at 
the  target  in  the  presence  of  wind  is  possible,  the  termi¬ 
nal  guidance  algorithm  determines  the  final  accuracy 
Homing  and  energy  management  are  necessary,  but 
not  critical,  in  determining  final  accuracy  because  it 
is  relatively  easy  to  bring  the  system  to  the  target,  but 
quite  difficult  to  plan  terminal  guidance  to  co-ordinate 
arriving  at  the  target  and  ground  simultaneously 

Examples  of  simple  strategies  are  Onyx  by  Atair 
Aerospace  [3,  4]  and  Screamer  by  Strong  Enterprises 
[3]  where  homing  takes  the  system  to  the  target, 
energy  management  is  a  spiral  around  the  target  at 
a  varying  radius,  followed  by  terminal  guidance  as  a 
straight  path  near  the  ground.  In  theory,  this  simple 
strategy  limits  the  error  to  the  spiral  radius.  How¬ 
ever,  in  practice  the  strategy  is  limited  to  the  case 
where  wind  speed  is  significantly  less  than  airspeed  in 
order  to  remain  over  the  target  and  maintain  guidance 
stability  [5]  making  it  unsuitable  for  high  wind-to- 
airspeed  ratios.  Another  limitation  of  the  spiraling 
strategy  is  that  it  cannot  guarantee  landing  into  the 
wind,  which  is  preferred  because  it  minimizes  speed 
at  impact.  Several  other  types  of  guidance  concepts 
for  parafoils  have  also  been  developed,  including  opti¬ 
mal  control,  trajectory  databases,  model  predictive 
control,  and  direct  glide  slope  control  [6-10].  More 
recently,  terminal  guidance  strategies  using  optimal 
trajectories  have  been  investigated.  Carter  etal.  devel¬ 
oped  a  method  for  generating  bandwidth  limited 
optimal  terminal  guidance  trajectories  [11].  Limiting 
the  trajectory  bandwidth  results  in  trajectories  that 
are  easily  tracked.  The  trajectory  was  parameterized 
using  five  parameters  with  the  optimization  problem 
iteratively  solved  to  minimize  final  position  and  head¬ 
ing  error.  Rademacher  et  al.  [12]  proposed  a  hybrid 
approach  where  optimal  trajectories  were  found  using 
either  modified  Dubins  paths  or  minimum  control 
trajectories  with  multiple  optimization  parameters. 
The  hybrid  strategy  was  employed,  because  each 
method  had  different  initial  conditions  where  nume¬ 
rical  robustness  was  a  problem  or  non-existence  of 
solutions  existed.  In  both  references  [11]  and  [12], 
kinematic  parafoil  models  were  used  and  the  optimal 
trajectories  were  tracked  using  proportional-integral- 
derivative  (PID)  controllers. 

This  article  develops  an  optimal  terminal  guidance 
algorithm  that  is  robust  in  high  winds  even  if  the 
wind  speed  exceeds  the  parafoil’s  airspeed.  An  opti¬ 
mal  terminal  trajectory  is  found  using  a  kinematic 
model  similar  to  references  [11]  and  [12].  However, 
using  the  direct  method  of  calculus  of  variations  and 
inverse  dynamics  in  a  virtual  domain,  the  problem  is 
reduced  to  a  single  parameter  optimization  problem 


[13] .  The  result  is  an  extremely  efficient  solution  that  is 
easily  solved  in  real-time.  In  addition,  since  the  entire 
terminal  trajectory  is  solved,  a  model  predictive  con¬ 
troller  (MPC)  [8]  is  used  to  track  the  trajectory  rather 
than  a  PID  algorithm.  A  complete  six  degrees-of- 
freedom  (DoF)  model  is  used  to  evaluate  the  method 
in  simulation  and  demonstrate  its  application  includ¬ 
ing  parafoil  turning  dynamics.  Finally,  experimental 
results  are  shown  where  the  algorithm  was  used  on  a 
small  2.3  kg  parafoil. 


2  GUIDANCE  STRATEGY 

In  order  to  ensure  the  parafoil  system  can  reach  a 
desired  target  it  is  released  upwind  of  a  target  at  high 
altitudes.  The  result  is  typical  parafoil  precision  place¬ 
ment  algorithms  having  an  energy  management  and 
homing  region.  The  proposed  strategy  for  precision 
placement  shown  in  Fig.  1  is  divided  into  four  phases 
with  the  first  being  energy  management,  the  second 
homing,  followed  by  terminal  guidance  comprising  a 
final  turn  and  final  approach  (FA).  Guidance  geometry 
is  defined  by  the  target  location  and  four  parameters: 
an  away  distance,  cycle  distance,  turn  diameter,  and 
wind  heading  angle. 

Energy  management  is  commenced  after  the  sys¬ 
tem  is  released.  The  parafoil  travels  clockwise  around 
the  rectangular  region  defined  by  the  away  distance, 
cycle  distance,  and  turn  diameter.  The  away  distance 
is  selected  so  that  in  the  case  of  wind  speed  exceed¬ 
ing  airspeed  the  parafoil  will  still  be  able  to  reach 
the  target,  while  in  energy  management,  the  parafoil 
continually  estimates  descent  rate  and  wind  speed. 
Descent  rate  and  wind  speed  estimates,  combined 
with  a  distance  from  the  target,  are  used  to  determine 
an  altitude  zstart  where  energy  management  is  termi¬ 
nated  and  homing  commences.  In  the  homing  phase 


Energy 

Management  . 

Homing  - 

Final  Turn  - 

Final  Approach  — 

Target  - 


Away 

Distance 


Cycle 

Distance 


\ 

\ 

l 


•  Diameter 
Release  X 


Wind 


Fig.  1  Guidance  strategy 
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the  system  travels  toward  the  target  while  continuing 
to  estimate  wind  and  vertical  velocity.  All  estimated 
data  are  used  to  determine  a  turn  initiation  point  (TIP) 
at  a  distance  Awitch  downwind  of  the  target  where  ter¬ 
minal  guidance  will  begin.  Terminal  guidance  starts 
with  a  final  turn  followed  by  an  FA  into  the  wind. 


3  TERMINAL  GUIDANCE 

Development  of  energy  management  and  homing 
strategies  is  flexible  with  the  only  requirement  being 
that  the  parafoil  is  brought  near  the  TIP  with  sufficient 
altitude.  In  contrast,  terminal  guidance  is  the  most 
critical  stage  of  parafoil  precision  placement  with  a 
very  strict  time  limitation.  The  parafoil  can  be  slightly 
late  or  earlier  departing  to  and  arriving  at  all  other 
phases,  but  terminal  guidance  ends  sharply  at  landing 
requiring  special  precautions  to  be  made  in  building 
a  control  algorithm.  An  ideal  terminal  guidance  tra¬ 
jectory  is  described  in  Fig.  2  where,  iI}  and  jIy  are  axes 
of  the  standard  North-East-Down  inertial  reference 
frame,  and  iT>  and  jT  are  target  axes  aligned  with  the 
wind  W.  The  parafoil  is  a  distance  L  upwind  of  the 
target,  tstan  is  the  current  time,  and  R  is  the  final  turn 
radius. 

The  final  turn  is  initiated  at  the  TIP,  occurring  at 
time  t0,  a  distance  Awitch  past  the  target.  After  the 
TIP  is  reached  the  final  turn  is  defined  by  a  com¬ 
manded  yaw  (t)  occurring  from  time  to  until  the  FA  is 
reached  at  time  A  The  resulting  final  turn  time  is  then 
Am  =  h  —  t0 .  The  desired  FA  time  rapp  is  t2  -  h  where 
t2  is  at  touchdown.  Transition  from  the  final  turn  to  the 
FA  occurs  at  Ak¬ 
in  terminal  guidance  described  above,  the  target 
location  R}  and  a  desired  Tapp  are  specified,  while  the 
wind  W  is  estimated  online  and  the  altitude  z  and  L  are 
measured.  The  terminal  guidance  problem  can  then 
be  summarized  as  follows.  For  a  parafoil  in  the  pres¬ 
ence  of  wind  Wy  at  altitude  z,  and  a  distance  L  from 


Fig.  2  Terminal  guidance  manoeuver 


the  target,  find  the  distance  Awitch  to  the  final  TIP  for 
an  ideal  impact  at  t2. 

In  general,  the  dynamic  model  of  a  parafoil  is  com¬ 
plex  and  non-linear,  so  that  the  problem  can  only  be 
solved  numerically.  However,  some  assumptions  can 
be  made,  allowing  an  analytical  solution  to  be  used  as 
a  reference  trajectory.  These  assumptions  are  a  slow 
turn  rate,  so  that  the  roll  and  sideslip  angles  can  be 
ignored,  and  nearly  constant  descent  rate  Vv  and  hori¬ 
zontal  airspeed  Ph.  In  this  case  it  immediately  follows 
that 


-^start 


(1) 


The  problem  reduces  to  a  kinematic  model  repre¬ 
sented  by  three  components  of  the  ground  velocity  in 
the  target  axes 


X 

—W  +  Ph  cos  \ff 

y 

= 

Ph  sin  i// 

z 

where  x,  y ,  and  z  are  positions  in  the  target  frame.  Now, 
the  commanded  turn  ^(t)  can  be  chosen  as  any  func¬ 
tion  that  satisfies  the  boundary  conditions  \lr(t0)  =  0 
and  VKh)  =  — ^  with  respect  to  the  wind  while  satis¬ 
fying  the  constraint  that  the  distance  travelled  in  jT  is 
-2 R.  For  a  constant  turn  rate,  the  solution  to  the  final 
turn  time  Ttuvn  and  turn  rate  are 


7T  R 

V* 


(3) 

(4) 


Hence,  the  most  straightforward  algorithm  to  con¬ 
trol  the  final  turn,  is  to  track  the  turn  rate  -VhIR. 
Assuming  that  the  wind  W  is  constant  and  using 
the  commanded  turn  rate  ijsc  =  —  VhIRf  integration  of 
inertial  velocities  along  iT  and  jT  from  fstart  to  t2  yield 


rh 


rtz 


switch  T 


x  dt  + 


=  A 


'switch 


WTt 


turn 


-  V 


(L_ t  -Awitch  \ 


x  At 

-  (Vh  +  uorapp  =  o 

+  Pv  Am  +  Pv  Tapp  =  0 


(5) 

(6) 


Equations  (5)  and  (6)  can  then  be  solved  for  Awitch 
and  Tapp  resulting  in 


/V.2  -  W2\ 

Awitch  =  WTturn  +  (  - -  ) 

( —z  rF  L  +  WTtmn  ^ 
X  \  Vv  ~  tUm  _  Vh  ~  W  ) 


L  +  WTturn 
ZVh 


(7) 

(8) 


From  equations  (7)  and  (8)  it  can  be  seen  that  the 
higher  the  altitude  z,  the  larger  Awitch  and  Tapp  become. 
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Fig.  3  Simulated  guidance  in  ideal  conditions  (•  TIP  and  Fig.  4  Final  turn  trajectories  with  initial  condition  errors 
FA  W  =  —3.4  m/s,  ■  TIP  and  FAW  =  -7.7  m/s) 


As  the  parafoil  loiters  upwind  of  the  target,  zstart  the  alti¬ 
tude  at  which  to  turn  towards  the  target  can  be  found 
using  a  desired  FA  time  rdpp.  The  switching  altitude  to 
exit  energy  management  which  achieves  rdpp  is  then 
given  by  solving  equation  (8)  for  z,  resulting  in 


^start  —  (  Aim  T 


L  +  WTturn 
Vh-w 


+ 


^des\ 

vh-w  appy 


(9) 


Once  the  system  is  homing  towards  the  target  the 
goal  is  to  bring  the  system  to  the  point  defined  by  x  = 
Awitch  and  y  =  —2R  with  Dswitch  continually  estimated 
during  homing  using  equation  (7)  and  current  esti¬ 
mates  of  Vh,  Vy,  and  W.  A  desired  FA  time  rdpp  is  used 
in  equation  (9)  in  making  the  decision  to  exit  energy 
management.  However,  disturbances,  measurement 
error,  and  tracking  error  while  homing  will  alter  the 
ideal  terminal  guidance  resulting  in  an  actual  rapp 
being  estimated  during  homing  using  equation  (7). 

Figure  3  demonstrates  trajectories  starting  150  m 
upwind  of  the  target  when  there  are  no  disturbances 
and  the  commanded  yaw  rate  is  tracked  precisely 
for  a  horizontal  airspeed  of  6.82  m/s,  descent  rate  of 
3.05  m/s,  turn  radius  of  37.5  m,  and  a  desired  FA  time 
rdpp  of  7.5  s.  Trajectories  for  two  wind  magnitudes  are 
shown,  -3.4  and  7.7  m/s,  resulting  in  the  wind  speed 
being  50  per  cent  and  110  per  cent  of  the  airspeed. 
The  manoeuvre  in  winds  of  50  per  cent  of  the  air¬ 
speed  starts  at  an  altitude  of  110  m  resulting  in  the 
TIP  at  a  distance  Awitch  of  -33.3  m  and  the  FA  begin¬ 
ning  25.7  m  downwind  of  the  target.  In  contrast,  when 
the  wind  speed  is  110  per  cent  of  the  airspeed  the 
manoeuver  starts  at  an  altitude  of  78  m  with  the  TIP 
at  a  distance  Awitch  of  -137  m  and  the  FA  beginning 
6.0  m  upwind  of  the  target,  with  the  parafoil  travelling 
backwards  at  impact. 


In  practice,  sensor  errors,  uncertain  winds,  and 
imperfect  control  will  disturb  the  ideal  touchdown 
depicted  above.  Figure  4  shows  Monte  Carlo  simula¬ 
tions  with  W  =  —3.4  m/s,  where  it  is  assumed  error 
in  the  initial  horizontal  position,  heading,  and  wind 
occur  at  the  beginning  of  terminal  guidance.  The 
standard  deviation  was  assumed  to  be  6  m  in  each 
co-ordinate,  10°  in  heading,  and  0.5  m/s  in  wind  speed. 
As  seen  from  Fig.  4,  these  errors  result  in  spreading 
the  arrival  to  the  FA  point.  Therefore,  while  the  sim¬ 
plified  algorithm  may  provide  a  reasonable  strategy 
for  exiting  the  energy  management  phase,  an  effec¬ 
tive  terminal  guidance  strategy  must  compensate  for 
sensor  errors  and  varying  winds  during  the  final  turn. 


4  OPTIMAL  TERMINAL  GUIDANCE 

To  overcome  tracking  errors  and  unaccounted  winds, 
the  following  two-point  boundary- value  problem 
(TPBVP)  is  formulated.  Staring  at  t  =  0  with  the  state 
vector  defined  as  x0  =  [Xo,yo,^o]T>  the  system  influ¬ 
enced  by  the  last  known  wind  needs  to  be  brought 
to  another  point,  Xf  =  [(14  +  M0Tdpp,  0,  — tt]t  at  t  =  tf. 
Hence,  they  need  to  find  the  trajectory  that  satisfies 
these  boundary  conditions  along  with  a  constraint 
imposed  on  yaw  rate,  \\j/\  ^  xj/max>  while  finishing  the 
manoeuver  in  exactly  Tturn  seconds.  The  optimal  x\r  (t)  is 
then  tracked  by  the  guidance  unit  using  any  appropri¬ 
ate  control  algorithm.  Errors  and  unaccounted  winds 
will  not  allow  exact  tracking  of  the  calculated  optimal 
trajectory;  therefore,  computation  of  an  optimal  tra¬ 
jectory  will  be  updated  during  the  final  turn,  each  time 
starting  from  the  current  state,  requiring  the  system 
to  arrive  at  Xf  =  [(14  +  V^)Tdpp,  0,  -tt]t  in  Am  seconds 
from  the  beginning  of  the  turn.  The  remaining  time 
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until  arrival  at  the  FA  is  computed  as 

=  (10) 

where  T^1  is  the  final  calculated  approach  time  from 
equation  (8)  before  entering  the  final  turn. 

Development  of  the  optimal  trajectory  algorithm 
begins  by  recalling  from  equation  (2)  the  target  frame 
horizontal  trajectory  kinematics 


"x" 

-(14  +  W)T«*f 

y. 

r=rf 

0 

X 

1 

1 

X 

"0" 

y. 

r=rf 

0 

} 

y. 

r=rf 

_0_ 

While  final  condition  equation  (17)  will  be  constant, 
initial  conditions  will  reflect  the  current  state  of  the 
system  at  each  cycle  of  optimization.  Differentiation 
of  equation  (15)  two  times  with  respect  to  x  results  in 


X 

~—W  +  14  cos  f 

y. 

14  sin  ijr 

(ID 


Recognizing  that  if  the  final  turn  horizontal  trajec¬ 
tory  is  given  by  equation  (11)  the  yaw  angle  along 
this  trajectory  is  related  to  the  change  of  inertial 
co-ordinates  as 


=  tan  1 


y 

x  +  W 


(12) 


Differentiating  equation  (12)  provides  the  yaw  rate 
control 


y(x  +  W)  -  xy 
(A  +  W02+y2 


(13) 


required  to  follow  the  reference  final  turn  trajectory 
in  presence  of  constant  wind  W.  The  inertial  speed 
along  the  trajectory  will  also  depend  on  the  current 
yaw  angle 


VG  =  v^+y2  =  yv^  +  W2-  2VbWcosf  (14) 


Now,  following  the  general  idea  of  direct  methods 
of  calculus  of  variations  [13]  the  solution  of  the 
TPBVP  will  be  represented  as  functions  of  some  scaled 
argument  f  =  x/xf  e  [0;  1] 


x(f)  =  P\(t)  =  a\  +  a\  x  +  a\  x2  +  a\  r3 
+  b\  sin(itr)  +  b\  sm(2nx) 
y(r)  =  P2(r)  =  &q  H-  a\  f  +  a\  x2  +  a\  r3 
+  b\  sin(Ttf)  +  b\  sin(2Kr) 

(15) 

Coefficients  a]  and  b]  {rj  =  1,2)  are  defined  by  the 
boundary  conditions  up  to  the  second- order  deriva¬ 
tive  at  x  =  0  and  r  =  rf  (f  =  1).  According  to  the  prob¬ 
lem  formulation  and  equation  (11),  these  boundary 
conditions  are 


ry\P' (f)  =  a\  +  2 a%r  +  3al f 2  +  tt b\  cos(irf) 

-1-  2kZ?2  cos(2kt) 

r^P''(r)  =  2a\  +  §a\x  —  tt 2b\  sin(jrf) 

—  (2n)2bl  sin(2Ttr)  (18) 


Equating  equation  (18)  at  the  terminal  point  to 
the  known  boundary  conditions  equations  (16)  and 
(17)  yields  a  system  of  linear  algebraic  equations  for 
coefficients  a]  and  (, b]rj  =  1,2).  For  instance,  in  the 
x-co-ordinate 


ao 

x0 

1  0  0  0  0  0 

11110  0 

a\ 

Xf 

0  10  0  tt  2tt 

a\ 

X'o  tf 
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X'Tf 
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i 
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Resolving  equation  (19)  yields 


a\  =  xq,  a\  =  -(xq-  xf)  + 


(2x'q  +  Xf)x2 


1  _  Tf 

2  o  > 


(Xq  -  X")Tf2 
6 

.//  |  v.//\_2 


1  2(x0  Xf)xf  +  (Xq  +  Xf)xf 

b  = - - - - - 

4tc 

12(Xo  —  xf)  +  6(Xq  +  xi)xf  +  (Xq  —  x"o)x2 
m  = - - - - - 

2  24tt 


(20) 


The  only  problem  is  that  derivatives  in  equation  (20) 
are  taken  in  the  virtual  domain,  while  the  actual 
boundary  conditions  are  given  in  the  physical  domain. 
Mapping  between  the  virtual  domain  [0;  rf]  and  phys¬ 
ical  domain  [0;  tf]  is  addressed  by  introducing  a  speed 
factor  X 


X 

"x0‘ 

"x" 

~—W  +  Vh  cos  iJ/0~ 

y_ 

r=0 

yo_ 

} 

y. 

r=0 

14  sin  ip0 

ix 

~-foVh  sin  i /f o" 

y_ 

r=0 

_  i/0vh  cos  i/Q  _ 

(21) 


Using  this  speed  factor  the  authors  may  now 
compute  corresponding  derivatives  in  the  virtual 
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domain  using  the  differentiation  rules  valid  for  any 
time -variant  parameter  § 

k  =  $  =  W&  +  *$")  (22) 

Inverting  equation  (22)  yields 

£'  =  A r1!,  §"  =  X~2%  -  X'X~^  (23) 

needed  to  transfer  the  boundary  conditions.  The  speed 
factor  X  simply  scales  the  entire  problem  -  the  higher 
speed  factor  1  is,  the  larger  zf  it  results  in  reference 
[14]  -  one  may  let  X  =  1  and  X'  =  0  at  the  boundary 
conditions  which  means  one  can  safely  assume  at 
the  boundary  conditions  £'  =  §  and  £ "  = 

Finding  the  optimal  solution  among  all  candidate 
trajectories  described  by  equation  (15)  begins  by  first 
guessing  a  value  of  the  only  varied  parameter  rf.  For 
the  initial  value  of  tf  one  can  take  the  length  of  the 
circumference  connecting  terminal  points  as 

rf  =  \^xj  -  x0)2  +  (y/  -  y0)2  (24) 


Time  (s) 

Fig.  5  Comparison  of  optimal  and  constant  turn  rate 

When  all  parameters  are  computed  at  each  of  the 
N  points,  one  can  compute  the  performance  index 

/  =  ^  Ar;  -  rturnj  +k*  a  (3i) 

where 

A  =  max  (0;  |t/^|  max)  (32) 


Coefficients  of  the  candidate  trajectory  are  com¬ 
puted  using  equation  (20)  with  the  boundary  condi¬ 
tions  equations  (16)  and  (17)  converted  to  the  virtual 
domain.  Having  an  analytical  representation  of  the 
candidate  trajectory  equations  (15)  and  (18)  defines 
the  values  of  xj}  yj}  x'j}  and  yj,  j  =  1, . . . ,  N  over  a  fixed 
set  of  N  points  spaced  evenly  along  the  virtual  arc  [0;  rf  ] 
with  the  interval 


At  =  Zf(N  —  1)  1 


(25) 


so  that 

Ti  =  ri- 1  +  At>  j  =  2, ...  ,N,  (n  =  0) 
Next,  for  each  node  j  =  2, . . . ,  N  compute 


At/_ i  = 


(xj-xj-i)2  +  (y]-y]-i)2 

K  +  W2  -  2V™  cos  fh  i 


(^i  =  to),  and 


(26) 


(27) 


a j  =  ArA  tf_\  (28) 

The  yaw  angle  f  can  now  be  computed  using  the 
virtual  domain  version  of  equation  (12) 


x/fj  =  tan 


-l 


XjXj  T  W 


Finally,  the  yaw  rate  \j/  is  evaluated  simply  as 

tj  =  (tj  - 


(29) 


(30) 


with  k+  being  a  weighting  coefficient. 

The  optimal  trajectory  problem  based  on  the  direct 
method  has  now  been  turned  into  a  single-parameter 
optimization  problem  where  zf  is  found  such  that 
equation  (31)  is  minimized  subject  to  the  system  in 
equations  (2)  and  (11).  The  problem  is  solved  using 
a  non-gradient  optimization  numerical  algorithm 
based  on  the  golden  section  search  and  parabolic 
interpolation  method  [15].  Reduction  of  finding  an 
optimal  terminal  trajectory  to  a  single-parameter  opti¬ 
mization  results  in  short  computation  times.  To  be 
more  specific,  a  1 6-bit  80  MHz  microprocessor  allowed 
computation  of  a  17.5  s  turn  manoeuver  with  N  =  25 
in  ten  iterations,  in  only  0.07  s. 

As  a  demonstration  of  the  proposed  approach,  Fig.  5 
shows  an  example  of  the  optimal  final  turn  for  the  case 
in  Fig.  3  when  W  =  -3.4  m/s.  The  optimal  final  turn 
and  constant  turn  rate  commands  are  similar;  how¬ 
ever,  as  opposed  to  the  constant  turn  rate,  the  yaw 
angle  changes  smoothly  at  departure  from  the  TIP  and 
arrival  at  the  FA.  As  a  result,  the  optimal  turn  trajectory 
is  slightly  different  but,  more  importantly,  the  system 
still  captures  the  FA  in  exactly  rtUrn  seconds.  Therefore, 
the  authors  have  a  tool  allowing  construction  of  the 
optimal  trajectory  from  any  initial  point  to  the  pre¬ 
determined  final.  To  this  end,  Fig.  6  shows  that  using 
the  same  initial  errors  from  Fig.  4  the  optimal  turn  is 
found  so  that  the  systems  still  arrive  at  the  FA  in  the 
predetermined  time,  rturn. 

Implementation  of  the  algorithm  on  actual  parafoil 
systems  requires  three  additional  components  to 
make  it  more  robust  in  practice:  compensation  for 
full  dynamic  motion,  wind  disturbance  mitigation, 
and  terminal  trajectory  tracking  errors.  The  kinematic 
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Fig.  6  Optimal  guidance  with  an  error  in  the  TIP 


model  in  equation  (2)  does  not  account  for  parafoil 
turning  dynamics  and  assumes  that  the  sideslip  and 
roll  angles  are  small.  When  the  turn  rate  is  sufficiently 
small,  or  the  radius  R  is  large,  the  model  in  equation  (2) 
provides  sufficient  accuracy  However,  in  order  to  track 
the  desired  trajectory  for  a  wide  range  of  R}  the  error 
from  sideslip  and  turning  dynamics  can  be  compen¬ 
sated  by  adding  an  additional  commanded  0C  for 
the  first  rpre  seconds  after  the  TIP  has  been  reached. 
Specifically,  for  the  first  fpre  seconds  of  the  maneuver 
the  command  yaw  0C  produced  by  equation  (29)  is 
augmented  as 

r  =  03) 

where  Kturn  is  the  correction  gain  and  Vh/R  is  the 
required  constant  turn  rate  in  equation  (4).  Distur¬ 
bances  from  wind  and  sensor  errors  during  tracking 
of  the  optimal  trajectory  will  result  in  errors  in  the 
FA  arrival.  These  errors  can  be  accommodated  by 
updating  the  optimal  trajectory  from  the  current  con¬ 
ditions  to  the  desired  FA  point  during  the  final  turn. 
Re-computation  of  the  optimal  trajectory  is  scheduled 
so  that  a  specified  number  of  updates  occur  at  equal 
intervals  over  the  period  from  t0  to  texit.  Finally,  the 
FA  location,  Xf  in  equation  (17),  assumes  that  during 
the  FA  the  parafoil  travels  directly  to  the  target  by  the 
shortest  path.  In  practice,  disturbances  and  the  guid¬ 
ance  algorithm  results  in  a  trajectory  with  less  than 
perfect  efficiency  at  reaching  the  target.  In  order  to 
correct  for  realistic  errors  during  the  FA  the  location  Xf 
is  corrected  to  be 

Xf  =  VhT«™  fiflnai  (34) 

where  £finai  is  the  FA  efficiency. 


5  PRECISION  PLACEMENT  SIMULATION 

In  the  previous  development  of  an  optimal  terminal 
trajectory  a  kinematic  model  was  used  and  it  was 
assumed  the  trajectory  was  tracked  perfectly.  How¬ 
ever,  actual  parafoil  systems  have  turning  dynamics 
resulting  in  side  slip  and  varying  descent  rates.  Evalu¬ 
ation  of  the  proposed  algorithm  will  be  done  using  a 
higher  fidelity  model.  The  combined  parafoil  canopy 
and  payload  can  be  modelled  using  six  DoF,  which 
include  three  inertial  position  components  of  the  sys¬ 
tem  mass  centre  as  well  as  the  three  Euler  orientation 
angles  of  the  parafoil-payload  system,  yaw  0,  pitch  0 , 
and  roll  cp.  The  6  DoF  model  described  in  reference  [9] 
is  used  in  all  simulations  with  a  constant  canopy  inci¬ 
dence  angle  r,  payload  drag  included  in  the  system 
aerodynamic  coefficients,  and  the  canopy  apparent 
mass  approximated  by  assuming  the  canopy  has  two 
axes  of  symmetry. 

The  guidance  strategy  outlined  in  sections  2  to  4 
requires  a  controller  that  tracks  a  commanded  yaw  0 c, 
using  the  brake  deflection  8a.  A  simple  two  DoF  model 
of  the  roll  and  yaw  dynamics  can  be  used  for  develop¬ 
ing  a  yaw  controller  since  for  parafoils,  pitch  and  speed 
are  not  typically  controllable.  The  state  vector  for  the 
two  DoF  rotational  model  is  x  =  [0,  0,p,  r]T,  where 
p  and  r  are  the  parafoil  roll  and  pitch  rates.  A  wide 
range  of  trajectory  tracking  controllers  used  on  UAVs 
could  be  used  [8,  16],  however,  MPC  is  selected  since 
the  optimal  terminal  guidance  prescribes  a  desired  0C 
over  the  terminal  trajectory  horizon.  The  rotational 
kinematics  and  dynamics  described  in  reference  [9] 
are  non-linear  and  thus  linearized  to  take  advantage  of 
well- developed  control  techniques  for  linear  systems. 
Assuming  that  the  aerodynamic  velocity  14  is  con¬ 
stant,  the  non-linear  model  with  parameters  listed  in 
Appendix  2  can  be  numerically  linearized  resulting  in 
a  linear  discrete  model  for  the  parafoil  with  sampling 
period  of  0.5  s 


<T 

0.962  0  0.153  0.012 

\jf 

0.0078  1  -0.011  0.043 

0 

P 

-0.103  0  0.033  0.004 

V 

r 

k+ 1 

0.0191  0  -0.0023  -0.003 

r 

-0.006 


_  0.1098  _ 

For  a  given  desired  discrete  output  vector  fc  = 
[0^+i,  0 k+2 >  •  •  •>  0 k+H  ]T  sPans  the  prediction  hori¬ 

zon  Hp)  MPC  solves  for  the  discrete  optimal  control 
vector  U  =  \uk,uk+i,...,uk+Hp-i]T  governed  by  the 
linear  plant  xk+i  =  Axk  +  B  uk  with  output  yk  =  Cxk 
[8,  17].  The  single-input  single-output  parafoil  model 
above  has  the  output  matrix  C  =  [010  0].  The  discrete 
linear  model  can  be  used  to  estimate  the  future  state 
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of  the  system  as 


—  KCa  Xk  U 


(36) 


from  the  current  state  vector  xk)  and  the  control  vector 
U}  where 


Kca 


Kcab  — 


CA 

CA2 

CAHp 


CB 

CAB 

CA2B 


CAHplB 


0 

CB 

CAB 


0 

0 

CB 


0 

0 

0 


0 

0 

0 

0 


CA2B  CAB  CB 


(37) 


(38) 


In  MPC,  the  optimization  problem  is  cast  as  a  finite¬ 
time  discrete  optimal  control  problem.  To  compute 
the  control  input  at  a  given  time  instant,  a  quadratic 
cost  function  is  minimized  through  the  selection  of 
the  control  history  over  the  control  horizon.  The  cost 
function  can  be  written  as 

/  =  (V  -  f  :)T  Q  (fc  -  f :)  +  UJRU  (39) 

where  Q  and  R  are  symmetric  positive  semi-definite 
matrices  of  size  Hp  x  Hp  that  penalize  tracking  error 
and  control,  respectively.  The  control  U}  which  mini¬ 
mizes  equation  (39),  can  be  found  analytically  [17]  as 

U  =  (K^abQKcab  +  R)  1  K£abQ  (r  -  KCAxk)  (40) 

Equation  (40)  contains  the  optimal  control  inputs 
over  the  entire  control  horizon.  In  practice,  only  the 
first  control  is  applied  and  the  problem  is  updated  at 
the  next  time  update. 

Implementation  of  the  terminal  guidance  algorithm 
requires  both  determining  the  optimal  terminal  tra¬ 
jectory  using  the  method  outlined  in  section  4  and 
the  trajectory  tracking  controller  in  equation  (40).  The 
first  step  in  the  complete  terminal  guidance  algorithm 
is  to  use  the  current  measured  parafoil  position  and 
velocity  along  with  the  known  FA  point  to  define  the 
boundary  conditions  to  the  TPBVP  (equations  (16), 
(17)  and  (34)).  Once  the  TPBVP  has  been  defined,  the 
guidance  algorithm  must  determine  the  optimal  solu¬ 
tion  among  all  candidate  trajectories  equation  (15)  by 
initializing  the  only  varied  parameter  Zf  according  to 
equation  (24).  The  optimal  trajectory  is  found  using 
the  single-parameter  optimization  problem  outlined 
in  section  4  (equations  (24)  to  (32))  with  a  non¬ 
gradient  optimization  numerical  algorithm  based  on 
the  golden  section  search  and  parabolic  interpolation 
method  [15].  Upon  completion  of  the  optimization 
problem,  the  desired  yaw  angle  trajectory  is  known 


from  equation  (29).  Finally,  the  discrete  optimal  yaw 
angle  can  be  combined  with  the  correction  in  equation 
(33)  and  used  in  the  proposed  MPC  algorithm  in 
equation  (40)  to  determine  the  final  control. 


6  SIMULATION  RESULTS 

In  all  simulations  the  full  non-linear  parafoil  model 
is  numerically  integrated  using  a  fourth- order  Runge- 
Kutta  algorithm  with  time  step  of  0.05  s.  The  MPC 
control  algorithm  described  previously  is  used  to  track 
Figures  7  and  8  show  examples  of  the  complete  ter¬ 
minal  guidance  strategy  for  an  R  of  37.5  m  in  winds  of 
-3.4  and  -7.7  m/s.  In  both  cases  rapp  is  7.5  s,  Sfinai  is 
0.95,  ti  —  is  3.0  s,  fpre  is  6.0  s,  Ktuvn  is  1.0,  V^max  is 
20  deg/s,  and  Jc*  is  400,  with  two  trajectory  updates 
made  during  the  final  turn.  Results  are  shown  from 
the  TIP  (x0,yo,£o)>  which  are  (-33,  75,  and  -93  m)  in 


Fig.  7  Optimal  real- dynamics- augmented  guidance  for 
W  =  -3.7  m/s  (a  trajectory  update) 


Fig.  8  Optimal  real- dynamics- augmented  guidance  for 
W  =  -7.7  m/s  (A  trajectory  update) 
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Fig.  9  Monte  Carlo  simulation  dispersion 


Fig.  8  and  (- 137, 75,  and  -93  m)  in  Fig.  9.  In  both  cases 
the  parafoil  initially  tracks  the  optimal  trajectory;  how¬ 
ever,  errors  slowly  build  due  to  mismatch  in  the  actual 
dynamics  and  the  model.  At  each  update  a  new  tra¬ 
jectory  is  calculated  that  results  in  arrival  at  the  FA. 
Impact  errors  are  0.4  m  in  Fig.  7  and  0.5  m  in  Fig.  8. 

The  complete  precision  placement  algorithm  com¬ 
bines  energy  management,  an  exit  from  energy  man¬ 
agement  based  on  equation  (9),  estimation  of  a  TIP 
during  homing  using  equation  (7),  and  the  optimal 
turn  followed  by  the  FA.  Monte  Carlo  simulations  of 
100  drops  were  done  using  the  complete  precision 
placement  algorithm  with  an  away  distance  of  450  m, 
a  cycle  distance  of  125  m,  an  R  of  37.5  m,  and  the  tar¬ 
get  at  the  origin.  The  initial  position  and  altitude  were 
nominally  760  m  upwind  of  the  target,  at  an  altitude  of 
700  m,  both  normally  distributed  with  a  50  m  standard 
deviation.  Gaussian  noise  was  injected  into  measured 
position,  altitude,  and  inertial  sensors  as  summarized 
in  Table  1. 

In  addition  to  sensor  errors,  two  sources  of  wind 
variation  were  added  to  the  simulation,  varying  direc¬ 
tion  and  varying  ground  wind  magnitude.  Wind  is 
assumed  to  have  a  constant  magnitude  prior  to  ter¬ 
minal  guidance  but  varies  linearly  from  W  at  the  TIP 


Table  1  Error  statistics 


Parameter 

la 

GPS  bias  (m) 

2.0 

GPS  noise  (m) 

0.5 

Altitude  bias  (m) 

2.0 

Altitude  noise  (m) 

0.5 

Roll,  pitch,  and  yaw  bias  (degree) 

2.0 

Roll,  pitch,  and  yaw  noise  (degree) 

1.0 

u,  v,  and  w  bias  (m/s) 

0.1 

u,  v ,  and  w  noise  (m/s) 

0.2 

p,  q,  and  r  bias  (degree /s) 

1.0 

p,  q}  and  r  noise  (degree/ s) 

1.0 

to  W  +  WF  at  the  ground.  The  wind  W  is  normally 
distributed  with  4.75  m/s  mean  and  2.0  m/s  standard 
deviation  while  WF  is  1.5  m/s  (la).  For  the  example 
parafoil,  this  results  in  a  nominal  wind-to-airspeed 
ratio  of  0.70  with  25  per  cent  exceeding  1.0.  Prevailing 
wind  is  assumed  by  the  system  to  travel  down  range 
at  a  heading  of  zero  degrees  while  the  true  wind  head¬ 
ing  varies  15°(la)  in  direction.  Dispersion  results  are 
shown  in  Fig.  9.  The  resulting  circular  error  probable 
(CEP)  shown  by  the  circle  is  16.8  m  and  is  defined  by 
the  radius  which  includes  50  per  cent  of  the  impacts. 

7  EXPERIMENTAL  RESULTS 

Test  results  of  the  small  2.3  kg  parafoil  outlined  in 
Appendix  2  are  shown  in  Fig.  10  with  an  away  dis¬ 
tance  of  200  m,  cycle  distance  of  120  m,  R  of  30  m, 
and  the  target  at  the  origin.  The  parafoil  is  released 
330  m  down  range  and  -240  m  cross  range  at  an  alti¬ 
tude  of  610  m  above  ground.  Wind  was  from  310° 
with  magnitude  varying  from  4.9  to  6.7  m/s  resulting 
in  a  wind-to-airspeed  ratio  of  70  per  cent  to  98  per 
cent.  Energy  management  is  exited  at  21 1  m,  97  s  after 
release.  Homing  proceeds  until  the  parafoil  is  60  m 
upwind  of  the  target,  121s  after  release,  at  an  alti¬ 
tude  of  103  m.  The  final  turn  lasts  9.0  s  until  the  FA 
begins  58  m  above  ground.  The  impact  occurs  7.4  m 
away  from  the  target  146  s  after  release. 

8  CONCLUSIONS 

A  terminal  guidance  algorithm  was  developed  for 
parafoil  precision  placement  in  high  wind-to-airspeed 
ratios.  Using  the  direct  method  of  calculus  of  variations 
and  inverse  dynamics  in  a  virtual  domain,  the  TPBVP 
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is  reduced  to  single-parameter  optimization  problem. 
The  result  is  an  extremely  efficient  solution  that  is 
easily  solved  in  real-time.  It  was  shown,  that  assum¬ 
ing  a  constant  final  turn  rate,  resulted  in  an  analytic 
solution  that  could  be  used  for  calculating  guidance 
decisions  prior  to  the  final  turn.  Uncertainty  during 
terminal  guidance  from  wind,  sensors,  and  guidance 
resulted  in  the  system  arriving  only  approximately  at 
the  predicted  location.  Efficient  computation  allowed 
the  trajectory  to  be  recomputed  during  the  final  turn 
adding  robustness  to  changing  winds  and  tracking 
errors.  The  algorithm  was  verified  using  a  full  6  DoF 
parafoil  model  and  Monte  Carlo  simulations  result¬ 
ing  in  a  circular  error  probable  of  16.8  m.  The  final 
test  results  were  presented  for  a  small  2.3  kg  system 
in  strong  winds  where  the  final  error  was  7.4  m. 

©Authors  2011 
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APPENDIX  1 


Notation 

A,  B,  C 

A,B,C 

b 

c 

d 

-^switch 

HP 

iijp  ki 
t  Ict 

id 

-^turn 

L 

p,q,r 

P,  Q,R 

Q, R 


R 

SP 

^exit 

^pre 


tapp> 


'T'des 

app 


'r’final 

app 


discrete  system  state  space  matrices 

element  of  the  apparent  mass  matrix 

canopy  span 

canopy  mean  chord 

maximum  brake  deflection 

distance  of  the  turn  initial  point  with 

respect  to  the  target 

discrete  prediction  horizon 

inertial  frame  unit  vectors 

target  frame  unit  vectors 

turn  rate  penalty  coefficient 

final  turn  correction  gain 

distance  to  target  along  target  line 

parafoil  roll,  pitch,  and  yaw  rates  in  the 

body  frame 

elements  of  the  apparent  inertia  matrix 

positive  semi-definite  error  and  control 

penalty  matrices 

final  turn  radius 

parafoil  canopy  area 

exit  time  for  final  turn 

length  of  time  to  apply  final  turn 

correction 

final  approach  time  and  desired  final 
approach  time 

last  calculated  final  approach  time 
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rturn  final  turn  time 

uy  Vy  w  parafoil  velocities  expressed  in  the 

body  frame 

U  discrete  optimal  control  vector 

Va  canopy  airspeed 

VG  parafoil  ground  speed 

Vh>  Vv  parafoil  horizontal  and  vertical  speed 

Wy  Wv  target  frame  wind  speed  and  ground 

wind  speed 

x,  y ,  z  parafoil  inertial  positions  in  the  target 

frame 

Zstart  altitude  to  exit  energy  management 

T  canopy  incidence  angle 

<$a  asymmetric  brake  deflection 

£ final  final  approach  efficiency 

A  speed  factor 

r  virtual  time 

rf  virtual  time  at  completion  of  final  turn 

f  virtual  time  scaled  by  the  final  virtual 

time 

4>y0y\lf  parafoil  Euler  roll,  pitch,  and  yaw 

angles 

\/fc  commanded  final  turn  angle 

ijr max  maximum  desired  final  turn  rate 

fc}  f  desired  and  estimated  output  vectors 

APPENDIX  2 

Model  parameters 

System  mass  m  (kg)  2.3 

Steady-state  aerodynamic  6.82 

velocity  14  (m/s) 


Canopy  reference  area 
Sp(m2) 

Canopy  span  b  (m) 
Canopy  chord  c  (m) 
Incidence  angle  r 
(degree) 

Inertia  matrix  elements 
(kg  m2) 

Elements  of  the 
apparent  mass  matrix 
(kg) 

Elements  of  the 
apparent  inertia 
matrix  (kg  m2) 
x-distance  from  mass 
centre  to  apparent 
mass  centre  xBM  (m) 
z-distance  from  mass 
centre  to  apparent 
mass  centre  zBM  (m) 
Maximum  brake 
deflection  d  (m) 
Aerodynamic 
coefficients 


1.1 

1.35 

0.75 

-12 

Ixx  =  0.423,  Iyy  =  0.401, 
Izz  =  0.052,  Ixz  = 
0.027 

A  =  0.012,  B  =  0.032, 

C  =  0.423 

P  =  0.054,  Q  =  0.13, 

R  =  0.0024 

0.05 


-1.1 


0.13 

Cdo  =  0.25,  Cd« 2  =  0.12, 
Cyp  =  0.23,  Q,o  = 
0.091,  CLa  =  0.90 
Cmo  =  0.35,  Cma  = 
—0.72,  Cmq  =  —1.49 
Qp  =  -0.036,  Clp  = 
-0.84,  Qr  =  -0.082, 
CUa  =  -0.0035, 

Cnp  =  -0.0015,  Cnp  = 
-0.082,  Cnr  =  -0.27, 
Cnsa  =  0.0115. 
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